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We propose a new mechanism for formation of topological defects in a U(l) model with a local 
gauge symmetry. This mechanism leads to definite predictions, which are qualitatively different 
from those of the Kibble-Zurek mechanism of global theories. We confirm these predictions in 
numerical simulations, and they can also be tested in superconductor experiments. We believe that 
the mechanism generalizes to more complicated theories. 
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When a global symmetry is spontaneously broken in a 
phase transition, it is generally accepted that the forma- 
tion of topological defects is well described by the Kibble- 
Zurek (KZ) scenario |^J^. As the transition takes place 
in a finite time, the correlation length of the order pa- 
rameter cannot keep up with its equilibrium value, which 
diverges at the transition point. The maximum correla- 
tion length reached determines the average distance of 
the defects in the final state , and it can be estimated 
from the critical dynamics of the theory . This scenario 
has been tested in many numerical simulations 1^,^ and 
in experiments with '^He Q and ^He Although the 
experimental results for ^He are in disagreement with the 
theory, the general picture is believed to be correct. 

However, in theories where the relevant symmetry is 
a local gauge invariance, e.g. in superconductors or in 
cosmology |7], the validity of the KZ scenario has been 
questioned | E ] , although lattice simulations ||] have found 
compatible results. Evidence that the KZ scenario might 
not work in these cases has been provided by recent ex- 
periments on YBCO superconductors [|lO|. Instead of 
detecting individual vortices, the authors measured the 
total net flux through the whole system, and found it to 
be zero within the accuracy of the experimental setup, in 
contradiction with predictions of the KZ scenario. 

The purpose of this letter is to suggest a simple and in- 
tuitive picture for defect formation in theories with local 
gauge symmetries. This picture is quite different from the 
KZ scenario and leads to definite predictions, which we 
have confirmed in numerical simulations. Furthermore, 
they can also be tested in superconductor experiments. 

Let us first review the KZ picture for the global case 
and consider for simplicity a U(l) symmetric field the- 
ory in D spatial dimensions. In the broken phase, the 
vacuum manifold is topologically a circle and the topo- 
logical defects are therefore vortices with dimensionality 
D — 2. The Hamiltonian contains a gradient term | V0p, 
where (f> = vexp{i6) is the order parameter field. Deep in 
the broken phase, the gradient term implies that V6' « 
in equilibrium. Thus the phase angle is correlated at in- 
finitely long distances, but since it takes an infinitely long 
time for the system to achieve that, after the transition 
6 will be approximately constant only inside domains of 



size ^, given by the maximum correlation length reached 
during the transition, and its value will be uncorrelated 
between these domains [Q. Since vortices are charac- 
terized by a non-zero change of 9 around a closed loop, 
they are formed where three domains meet with a prob- 
ability that is independent of the size of the domains. 
Consequently, the final vortex number in the global case 
behaves as A'' ~ The value of ^ can be estimated 

from the critical slowing down of the dynamics during 
the transition 

If the symmetry is local, there is another, competing 
mechanism, which dominates if the transition is suffi- 
ciently slow. In the temporal gauge — 0, the Hamil- 
tonian for a relativistic gauged U(l) scalar field theory 
in three spatial dimensions is 
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where H = do(j) is the canonical momentum, E — — 
and B = V X A the electric and magnetic field strengths, 
and V{4)) the potential of the scalar field. The corre- 
sponding equations of motion are 



dicj, = &cj> - y (0), 

DqE = V xB + 2eIm(j)*D(j}, 
V -E^ 2elm(/)*ao0. 
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(Note that we use units with ^ fi = c = iiq = 1.) 
More generally, we will consider the analogous system in 
D spatial dimensions. Because the gradient has been re- 
placed by a covariant derivative D(f) = \7(f) + ieAcj), the 
energy is minimized in the broken phase if V0 ~ —eA. 
In the presence of a magnetic flux, this condition can- 
not be satisfied everywhere and frustrations, vortices, are 
formed. Although the magnetic flux is zero on the aver- 
age, the thermal fluctuations give it a non-zero variance. 
When the system enters the broken phase, it tries to re- 
arrange the field configuration to minimize the magnetic 
flux and the energy associated with it. Because of the 
finite time available, this is not possible for the fluctu- 
ations with the longest wavelengths and they freeze in 
their initial form, but at shorter distances the fluctua- 
tions of the magnetic flux are smoothed out. Therefore, 
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immediately after the transition, before the flux is local- 
ized into vortices, its configuration consists of domains 
inside which it is approximately uniform and which have 
some characteristic size ^. 

If we calculate the winding number around a curve C, 
which encircles one of the domains discussed above, it 
typically does not vanish, but instead 
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where $0 is the magnetic flux through the curve C. Be- 
cause the flux rearrangements at the transition were only 
able to change its distribution inside the domain, $(7 h^^s 
the same value it had in the symmetric phase. We can 
estimate it in the standard way by calculating the energy 
E{^c) associated with the flux. Inside the domains, the 
flux is uniform and therefore 



(4) 



Requiring E{^c) ~ T gives the typical value of the flux 
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Using Eqs. (^) and (^, we can now estimate the area 
density of vortices after the transition. Suppose we have 
in our D-dimensional space a surface of area A, then it 
is split into ^ A£^^^ domains in the transition. Each 
domain contains Nq « {e/2T:)^c vortices, and we will 
assume that iVg > 1. Then, the total number of vortices 
piercing the surface per unit area is 
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In particular, at D = 2, iV ~ 

In addition, vortices will also be formed by a variant 
of the KZ mechanism. Initially, these vortices behave 
like global vortices, but eventually a quantum of mag- 
netic flux is generated inside each of them, making them 
truly local vortices. However, the KZ mechanism is only 
important if the transition is very rapid or the tempera- 
ture very low, i.e. Nq < 1. Thus we will neglect it in the 
following. 

At D = 3, the vortices formed inside the domains form 
a network at distances longer than ^, and in this network 
most vortices will be in the form of closed loops, which 
will quickly shrink into a point and disappear. Therefore 
it is useful to consider a border-line case between D = 2 
and _D = 3, where one of the dimensions Lz is very short. 
As long as Lz < the vortices will wind around the short 
dimension rather than forming loops and therefore they 
will be stable. Now the domains have the form of a short 
cylinder, and the estimate (^) generalizes easily to this 
case, yielding 



N/A « ^t1/2l;1/2^'-i _ ^-1/2^ 
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The correlations between the vortices in the final state 
will also be different from those predicted by the KZ 
mechanism. Let us consider a system after the transi- 
tion and assume that there is a vortex with a positive 
winding at point x. In the KZ scenario, the distance to 
the nearest other vortex should be roughly If we cal- 
culate the winding number nc{r) of a circular loop C of 
radius r centered at x, it follows that at r < ^, the wind- 
ing number is close to one. However, at distances r>^, 
the phase angle is independent of whether there is a vor- 
tex inside at x or not, and therefore nc{r) = 0. Thus, in 
the KZ scenario, 
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Again, our case is very different. Inside a single do- 
main, all vortices have the same sign, and ncir) is there- 
fore an increasing function of r at r < ^. At r > ^, nc{r) 
gets contributions from different domains, and since they 
are all independent they average to zero and nc{r) be- 
comes a constant. The behaviour in our scenario is there- 
fore totally different from Eq. M): 
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where cq > and ci > I are constants. 

Finally, let us discuss the dependence of ^ on the 
"quench" timescale tq, which parameterizes the rate at 
which the phase transition takes place. For definiteness, 
we consider the potential 
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where the mass parameter is changed linearly across its 
critical value 
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In reality, is not equal to zero when thermal fluctua- 
tions are taken into account, but let us use the mean- field 
approximation in the following, in which case it is. 

If we know the photon dispersion relation lo — u}{k) 
at the transition point and its neighbourhood, we can 
estimate that a Fourier mode of wave number ki falls out 
of equilibrium during the transition if the adiabaticity 
condition is not satisfied 
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However, the calculation of the dispersion relation is 
beyond the scope of this letter, and we will instead 
only consider two simple special cases: the overdamped 
case (OD), where the dynamics is dominated by a k- 
independent damping rate 7, and the underdamped case 
(UD) with the free-field dispersion relation: 
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i(fc2+m2) (OD), (fc2 + m2)i/2 (UD). (13) 
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These same special cases have been discussed in the con- 
text of global theories in Ref. but we stress that 
it is by no means clear that the dynamics is well de- 
scribed by either of these cases. Furthermore, we assume 
that the photon mass for the relevant modes behaves as 
« 2e^|0p, and that is given by its equilibrium 
value |0p = -m^{t)/2\ - t/rg. 

Now, Eq. ( p^ ) tells us that the highest wave numbers 
that fall out of equilibrium behave as 
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The domain size ^ is then simply given by ^ w 27r/fc, and 
using Eq. (|^), we can write down the dependence of the 
final vortex number on tq in two dimensions as 
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In the KZ scenario, the analogous exponents are 
and —2/3, respectively 

We carried out a set of numerical simulations using the 
equations of motion ^ to test the results in Eqs. (|^), 
(^ and (|l^). Our coupling constants were e = 0.3 and 
A = 0.18. Since A > e^, the transition is continuous, as 
in a Type II superconductor. We used periodic lattices 
of size 120 X 120 x L^, where = 5 or 20 (in units 
where the lattice spacing is one). We prepared a set of 
initial conditions according to the thermal distribution 
exp{—H/T) at T = 6 using a hybrid Monte Carlo algo- 
rithm and followed the time evolution using a leap-frog 
algorithm with time step St = 0.05, changing the mass 
parameter according to 
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Here mg = —1.6 is the initial value at t = — tq when the 
quench begins and 5m? — 3.2 was chosen in such a way 
that the transition takes place at t « 0. As in Eq. (11), 
the rate of change of m? is proportional to tq , but in this 
form, we reach an equilibrium state at late times, which 
makes comparison of different quench rates easier. The 
details of the simulations will be discussed in Ref. . 

In the final state, the vortices have only short-range 
interactions and therefore they freeze quickly into their 
final configuration. After that, at t = tq -I- 400, we lo- 
cated the vortices from the field configuration by mea- 
suring the gauge- invariant winding numbers of individual 
plaquettes ]l3[ | and connecting plaquettes with non-zero 
winding numbers into vortex lines. We then measured 
the total number of those vortex lines that wind around 
the short dimension of our lattice. The results are shown 
in Fig. 1^ as a function of tq. The results from the thinner 
lattice = 5 have been divided by 2, because according 
to Eq. (0) the results from the two lattices should then 
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FIG. 1. The number of vortices in the final state as a func- 
tion of TQ. The Lz = 5 data have been scaled by a factor of 
1/2 in accordance with Eq. (^. The solid and dashed lines 
are the fits to Eq. (^) (see Table ^ and correspond to the 
cases = 20 and — 5, respectively. 
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Tc 


a 


x' 




= 5 


49.5 ±2.6 


8.5 ± 1.7 


0.250 ±0.013 


6.73 




= 20 


31.2 ±6.7 


54.1 ±13.7 


0.274 ±0.039 


12.09 



TABLE I. Fit of the data to the functional form (|l7|). Re- 
sults from both simulations are compatible with the over- 
damped exponent a = 0.25 (see Eq. dlq)). 



be on top of each other. Each datapoint is an average 
of ^ 20 runs starting from different initial configurations 
drawn from the thermal ensemble. 

In slow transitions, the data agree very well with the 
predicted power-law behaviour. At tq < 40— 50, the — 
20 data become independent of tq, which suggests that 
^ < Lz, and part of the vortices have formed loops and 
do not contribute to the vortex count. We fitted the 
function 
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into the data, and the results are given in Table |. Both 
cases agree with the overdamped prediction a — —0.25. 
Moreover, the ratio of the results from — 5 and — 
20 measured at, say, tq — 100, is 1.8 ± 0.2, which is 
compatible with the value 2 predicted by Eq. (|^). 

Using our data, we can also perform a more quanti- 
tative check for Eq. (^ even without knowledge of the 
dynamics that determines ^, if we assume that ^ = Lz 
at the turning point in the = 20 data. By substitut- 
ing this into Eq. (|^) , we find that the number of vortices 
formed at tq < Tc in the Lz — 20 case should he N fa 19, 
which agrees within a factor of 2 with the numerical re- 
sult. 

We can also test Eq. (^ by measuring the net num- 
ber of vortices inside a circle of radius r centered at a 
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' Xq=500, t=500 
. Xq=10, t=240 
random 




FIG. 2. The net number of vortices within distance r from 
a vortex with positive winding in a slow (circles) and a fast 
(triangles) quench on a thin lattice [Lz = 5). The dotted line 
is the corresponding curve for uncorrelated vortex-antivortex 
pairs. The fact that the data points are well above the ran- 
dom curve shows that there is a positive correlation between 
vortices at short distances. 



which is below the resolution of the experiment. A sim- 
ilar estimate applies for a recent Josephson junction ex- 
periment where iV « 7 flux quanta were observed. 
In this case, the prediction of the KZ scenario is TV « 4, 
and thus the experiment cannot decide between the two 
mechanisms. It seems that our scenario can only be con- 
firmed by experiments in which not only the total flux 
but also the spatial distribution of the vortices can be 
measured. 

To summarize, our numerical simulations confirm the 
results in Eqs. (|^) and (^, which are independent of as- 
sumptions about the dispersion relation u}{k). Using the 
naive overdamped dispersion relation (^3|), we can also 
reproduce accurately the exponents a in Table ||. This 
supports strongly the scenario presented in this letter. 
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vortex. The results for two different values of tq with 
Lz = 5 are shown in Fig. |^, together with the benchmark 
curve nc = 1 — nr'^A^^, which corresponds to a uni- 
form random distribution of vortcx-antivortex pairs and 
shows the effect of the finite system size. At short dis- 
tances, the data points are significantly above this curve, 
which indicates a positive correlation between vortices, 
and at long distances, they follow the benchmark curve. 
Both of these results agree with Eq. (||) but differ from 
the KZ prediction (|). 

Our results disagree with the simulations in Ref. 
where the authors found an exponent a that was compat- 
ible with the global theory. We believe that the temper- 
ature used there, T — 0.01, was so low that practically 
all of the vortices were formed by the KZ mechanism. 

Since our mechanism can only increase the number of 
vortices formed in a transition from that predicted by the 
KZ mechanism, it cannot explain the failure in Ref. ||ic| ] 
to find any total net magnetic flux when a superconduc- 
tor film was quenched into the superconducting phase. 
However, the experiment does not rule it out either, be- 
cause the vortices formed by the KZ mechanism might 
avoid detection by being expelled from the film before 
generating observable magnetic flux. The extra magnetic 
flux predicted by our mechanism is fairly small, because 
it can only change the flux distribution at short distances. 
The minimal energy for a configuration with a given value 
$ of fiux through the film is that of a magnetic dipole, 
Emin ~ n^^A-^/'^^'^. Using the values T = 90 K and 
A = 1 cm^, we find that the predicted number of flux 
quanta is 
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